---
title: "Lab 8 Code"
author: "Fabio Setti"
date: "Spring-2025"
output: html_document
---


## Today's Packages and Data 🤗


```{r eval=FALSE}
# set to eval = TRUE to install packages

# install.packages("tidyverse")
install.packages("lubridate")
install.packages("splines2")
```


```{r}
library(tidyverse)
theme_set(theme_classic(base_size = 14, 
                        base_family = 'serif'))
library(lubridate)
library(splines2)
```

```{r}
NY_waeth <- rio::import("https://github.com/quinix45/PSYC-7804-Regression-Lab-Slides/raw/refs/heads/main/Slides%20Files/Data/NYC_Weather_2016_2022.csv")

str(NY_waeth, vec.len = 2)
```

## Working with Dates


```{r}
NY_waeth$time[1:3]
```



```{r}
# substr() takes a sting and extract characters from (1 here) to (10 here)

NY_waeth$time <- substr(NY_waeth$time, 1, 10)
NY_waeth$time[1:3]
```


```{r}
NY_waeth$time <- lubridate::ymd(NY_waeth$time)
class(NY_waeth$time)
```


```{r}
NY_waeth$time[1]
as.numeric(NY_waeth$time[1])/365
```



## Summarizing and Filtering the Data


```{r}
Temp_avg <- NY_waeth %>% 
              select(time, `temperature_2m (°C)`) %>%
                group_by(time) %>% 
          summarise(Temp_C = mean(`temperature_2m (°C)`))

str(Temp_avg)
```


```{r}
Temp_avg$Temp_F <- (Temp_avg$Temp_C * 9/5) + 32 
```



```{r}

Temp_avg_2021 <- Temp_avg %>% 
                  filter(time >= "2021-01-01", 
                         time <= "2021-12-30")
str(Temp_avg_2021)
```



## Plot The Data


```{r}
ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point()
```



```{r}
ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point() +
geom_smooth(method = "lm", 
            formula = y ~ poly(x, 2), 
            se = FALSE)
```



```{r}
ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point()  +
geom_smooth(method = "lm", 
           formula = y ~ poly(x, 2), 
           se = FALSE) +
geom_smooth(method = "loess", 
            col = "red",
           se = FALSE)
```

## Quadratic Regression: Centering Variables


```{r}
Temp_avg_2021$time2 <- as.numeric(Temp_avg_2021$time)^2
```

```{r}
cor(Temp_avg_2021$time2, as.numeric(Temp_avg_2021$time))
```

```{r}
Temp_avg_2021$time_cent <- Temp_avg_2021$time - mean(Temp_avg_2021$time)
# btw, check out what the time_cent variable looks like now!

# we create the quadratic term by squaring the centered variable
Temp_avg_2021$time_cent2 <- as.numeric(Temp_avg_2021$time_cent)^2
```


```{r}
cor(Temp_avg_2021$time_cent2, 
    as.numeric(Temp_avg_2021$time_cent))
```



## Quadratic Regression


```{r}
quad_reg <- lm(Temp_F ~ time_cent + time_cent2 , 
               data = Temp_avg_2021)

summary(quad_reg)
```



## $b_0$ in Quadratic Regression


```{r}
ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point(alpha = .5, shape = 1)  +
  geom_point( aes(y = coef(quad_reg)[1], 
                  x = mean(Temp_avg_2021$time)),
              size = 3.5,
              col = "red") +
geom_smooth(method = "lm", 
           formula = y ~ poly(x, 2), 
           se = FALSE) +
geom_segment(x = 0, 
               xend = mean(Temp_avg_2021$time), 
               y = coef(quad_reg)[1],
               yend = coef(quad_reg)[1],
               linetype = 2) +
geom_segment(x = mean(Temp_avg_2021$time), 
             xend = mean(Temp_avg_2021$time), 
             y = 0,
             yend = coef(quad_reg)[1],
             linetype = 2)
```


## $b_1$ in Quadratic Regression


```{r}
ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point(alpha = .5, shape = 1)  +
geom_smooth(method = "lm", 
            formula = y ~ poly(x, 2), 
            se = FALSE) +
geom_segment(x = mean(Temp_avg_2021$time) - 130, 
             xend = mean(Temp_avg_2021$time) + 130, 
             y = coef(quad_reg)[1] - coef(quad_reg)[2]*130,
             yend = coef(quad_reg)[1] + coef(quad_reg)[2]*130,
             linetype = 1,
             linewidth = .7,
             col = "red")
```



## $b_2$ in Quadratic Regression




```{r}
x_val <- 100

# accurate calculations, there is rounding error in the equation on the slides
slope_new <- coef(quad_reg)[2] + 2*coef(quad_reg)[3]*x_val

ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point(alpha = .5, shape = 1)  +
geom_smooth(method = "lm", 
            formula = y ~ poly(x, 2), 
            se = FALSE) +
geom_segment(x = (mean(Temp_avg_2021$time) + x_val) - 130, 
             xend = (mean(Temp_avg_2021$time) + x_val) + 130, 
             y =  (coef(quad_reg)[1] + coef(quad_reg)[2]*x_val + coef(quad_reg)[3]*x_val^2) - (slope_new*130),
             yend = (coef(quad_reg)[1] + coef(quad_reg)[2]*x_val + coef(quad_reg)[3]*x_val^2) + (slope_new*130),
             linetype = 1,
             linewidth = .7,
             col = "red") 
```


## Climbing to the top (or bottom)


```{r}
# accurate calculations, there is rounding error in slide regression coefficients
-1*(coef(quad_reg)[2]/(2*coef(quad_reg)[3]))
```



```{r}
# accurate calculations, there is rounding error in the equation on the slides
top <- -1*(coef(quad_reg)[2]/(2*coef(quad_reg)[3]))

Y_val <- (coef(quad_reg)[1] + coef(quad_reg)[2]*top + coef(quad_reg)[3]*top^2)

# this is the formula for the tangent slope. if you print this, you will see that it equals exactly 0. This means that the X value provided is the top of the curve, which is the only point where the tangent slope is exactly 0

slope_new <- coef(quad_reg)[2] + 2*coef(quad_reg)[3]*top


ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point(alpha = .5, shape = 1)  +
  geom_point(aes(y = Y_val, 
                x = mean(Temp_avg_2021$time) + top),
              size = 3.5,
              col = "red") +
geom_smooth(method = "lm", 
           formula = y ~ poly(x, 2), 
           se = FALSE) +
geom_segment(x = 0, 
               xend = mean(Temp_avg_2021$time) + top, 
               y = Y_val,
               yend = Y_val,
               linetype = 2) +
geom_segment(x = mean(Temp_avg_2021$time) + top, 
             xend = mean(Temp_avg_2021$time) + top, 
             y = 0,
             yend = Y_val,
             linetype = 2) +
geom_segment(x = mean(Temp_avg_2021$time) -130, 
             xend = mean(Temp_avg_2021$time) + 130, 
             y =  Y_val,
             yend = Y_val,
             linetype = 1,
             linewidth = .7,
             col = "red") 
```


## Bootsrap CI for Maximum/Minumum point




```{r}
boot_quad <- car::Boot(quad_reg, R = 1000)
```


```{r}
# View(boot_quad$t)
# we can compute -b1/(2*b2) for all sample at once with this

top_boot <- -boot_quad$t[,2]/(2*boot_quad$t[,3])
```


```{r}
quantile(top_boot, c(.025, .975))
```


```{r}
plot(density(top_boot))
```


## How Many Polynomials do I need?


```{r}

reg_lin <- lm(Temp_F ~ time, data = Temp_avg_2021)  
reg_quad <- lm(Temp_F ~ poly(time, 2), data = Temp_avg_2021)  
reg_cube <- lm(Temp_F ~ poly(time, 3), data = Temp_avg_2021)  
reg_fourth <- lm(Temp_F ~ poly(time, 4), data = Temp_avg_2021) 
reg_fifth <- lm(Temp_F ~ poly(time, 5), data = Temp_avg_2021) 
```



## Comparing Polynomial Regressions


```{r}
c(summary(reg_lin)$r.squared,
  summary(reg_quad)$r.squared,
  summary(reg_cube)$r.squared,
  summary(reg_fourth)$r.squared,
  summary(reg_fifth)$r.squared)

```


```{r}
anova(reg_lin, reg_quad, reg_cube, reg_fourth, reg_fifth)
```


## Plotting polynomials

```{r}
ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point() +
geom_smooth(method = "lm", 
           formula = y ~ poly(x, 2), 
           se = FALSE) +
  ggtitle("Second Degree Polynomial") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```


```{r}
ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point() +
geom_smooth(method = "lm", 
           formula = y ~ poly(x, 3), 
           se = FALSE) +
  ggtitle("Third Degree Polynomial") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```



```{r}
ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point() +
geom_smooth(method = "lm", 
           formula = y ~ poly(x, 4), 
           se = FALSE) +
  ggtitle("Fourth Degree Polynomial") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```


```{r}
ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point() +
geom_smooth(method = "lm", 
           formula = y ~ poly(x, 5), 
           se = FALSE) +
  ggtitle("Fifth Degree Polynomial") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```



## Piecewise Regression



```{r}
# split data with cut-point (cleaner code in 2 slides, can't be bothered to change it)

Temp_avg_2021$D1 <- ifelse(Temp_avg_2021$time_cent >= 20.34353, 0*Temp_avg_2021$time_cent, 1*Temp_avg_2021$time_cent)
Temp_avg_2021$D2 <- ifelse(Temp_avg_2021$time_cent < 20.34353, 0*Temp_avg_2021$time_cent, 1*Temp_avg_2021$time_cent)

reg_piece <- lm(Temp_F ~ D1 + D2, data = Temp_avg_2021)

intercept <- coef(reg_piece)[1]
D1 <- coef(reg_piece)[2]
D2 <- coef(reg_piece)[3]


ggplot(Temp_avg_2021,
       aes(x = time, y = Temp_F)) +
        geom_point() +
geom_segment(x = min(Temp_avg_2021$time_cent) + mean(Temp_avg_2021$time), 
             xend = 20.34353 + mean(Temp_avg_2021$time), 
             y = intercept + min(Temp_avg_2021$time_cent)*D1 ,
             yend = intercept,
             linetype = 1,
             linewidth = .7,
             col = "red") +
geom_segment(xend = max(Temp_avg_2021$time_cent) + mean(Temp_avg_2021$time), 
             x = 20.34353 + mean(Temp_avg_2021$time), 
             yend = intercept + max(Temp_avg_2021$time_cent)*D2 ,
             y = intercept,
             linetype = 1,
             linewidth = .7,
             col = "red") +  
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))
```




## Setting up Piecewise 

```{r}
# an example of of ifelse() works

ifelse(c(2,3) > 2.5, "bigger", "smaller")
```



## Creating $D_1$ and $D_2$


```{r}
Temp_avg_2021$D1 <- ifelse(Temp_avg_2021$time_cent >= 20.34353, 
                           0, Temp_avg_2021$time_cent)
```



```{r}
Temp_avg_2021$D2 <- ifelse(Temp_avg_2021$time_cent < 20.34353, 
                           0, Temp_avg_2021$time_cent)
```

## Running The Piecewise Regression

```{r}
piece_reg <- lm(Temp_F ~ D1 + D2, data = Temp_avg_2021)
summary(piece_reg)
```


## 3D visualization? Out of Curiosity...


```{r}
library(FabioFun)

nice_3D_plot(y = Temp_avg_2021$Temp_F,
             x1 = Temp_avg_2021$D1,
             x2 = Temp_avg_2021$D2,
             axis_names = c("Temperature", 
                            "D1", 
                            "D2"),
             reg_plane = TRUE)
```



## Spline Regression


```{r}
spline_fit <- lm(Temp_F ~ splines2::bsp(Temp_avg_2021$time_cent, df = 5), 
                 data = Temp_avg_2021)
summary(spline_fit)
```


## Visulizing Spline Regression

```{r}

pred_data <- seq(min(Temp_avg_2021$time_cent),
                 max(Temp_avg_2021$time_cent))

pred <- predict(spline_fit, 
                newdata = list(time_cent = pred_data))

```



```{r}
ggplot(Temp_avg_2021,
       aes(x = time_cent, y = Temp_F)) +
        geom_point() +
  geom_smooth(method = "loess", 
              formula = y ~ x,
              se = FALSE) +
  geom_line(aes(x = pred_data,
                y = pred),
            col = "red",
            size = 1) 
```



